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We have implemented time-propagation of the non-equilibrium Green function for atoms and 
molecules, by solving the Kadanoff-Baym equations within a conserving self-energy approximation. 
We here demonstrate the usefulnes of time-propagation for calculating spectral functions and for 
describing the correlated electron dynamics in a non-perturbative electric field. We also demonstrate 
the use of time-propagation as a method for calculating charge-neutral excitation energies, equivalent 
to highly advanced solutions of the Bethe-Salpeter equation. 
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The arrival of molecular electronics has exposed the 
need for improved methods for first-principles calcula- 
tions on non-equilibrium quantum systems [l| . The non- 
equilibrium Green function 0, 0| is for several reasons a 
natural device in such studies. Not only is it relatively 
simple, being a function of two coordinates, but it con- 
tains a wealth of information, including the electron den- 
sity and current, the total energy, ionization potentials, 
and excitation energies. Time-propagation according to 
the Kadanoff-Baym equations [3, |4| is a direct method 
for describing the correlated electron dynamics, and a 
method which automatically leads to internally consis- 
tent and unambiguous results in agreement with macro- 
scopic conservation laws. In the linear response regime, 
time-propagation within relatively simple self-energy ap- 
proximations corresponds to solving the Bethe-Salpeter 
equation with highly advanced kernels The Green 
function techniques are also highly interesting as a com- 
plementary method to time-dependent density functional 
theory (TDDFT) 0,0], not only by providing benchmark 
results for testing new TDDFT functional, but diagram- 
matic techniques can also be used to systematically de- 
rive improved density functionals 0]. We will in this 
Letter demonstrate time-propagation for inhomogenous 
systems, using the beryllium atom and the H 2 molecule 
as illustrative examples. 

The non-equilibrium Green function G(xi, x'i') de- 
pends on two time- variables, rather than one such as 
the time-dependent many-particle wave function. On the 
other hand, the fact that it also depends on only two 
space and spin variables x = (r, a) means that it can 
be used for calculations on systems that are too large for 
solving the time-dependent Schrodinger equation, such as 
the homogenous electron gas [j| or solids Q , or for calcu- 
lations on molecular conduction. In addition, the Green 
function provides information about physical properties 
such as ionization potentials and spectral functions, that 
are not given by the many-particle wave function. The 
Green function techniques also have the advantage over, 
e. g., TDDFT or density matrix methods based on the 



FIG. 1: At t — 0, the system is in thermal equilibrium, and 
the Green function is calculated for imaginary times from 
to — i/3 (left figure). Describing the system for t > implies 
calculating G(tx,ti) on an extending time-contour (right). 



BBGKY hierarchy [13], that it is easy to find approxima- 
tions which give observables in agreement with macro- 
scopic conservation laws. 

The two time-arguments of the Green function are lo- 
cated on a time-contour as illustrated in Fig. [T] It solves 
the equation of motion (we suppress the space- and spin- 
variables for notational simplicity) 

[id t -h(t)]G(t,t') = S(t,t')+ { diT,(t,i)G(i,t'), (1) 
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where h(t) is the non-interacting part of the hamiltonian 
(which may include an arbitrarily strong time-dependent 
potential), and the self-energy £ accounts for the effects 
of the electron interaction. We use atomic units through- 
out. The time-integral is performed along the contour, 
and the delta- function 5(t,t') is defined on the contour 
[rH . We only consider systems initially (at t = Q) in the 
ground state. This is facilitated by describing the system 
within the finite-temperature formalism [T^, letting the 
time-contour start at t = and end at the imaginary time 
t = —ij3 — ~i/(ksT), as illustrated in Fig. [TJ This choice 
of initial conditions means that the first step consists of 
calculating the Green function for both time-arguments 
on the imaginary track. The calculations are carried out 
in a basis of Hartree-Fock (HF) molecular orbitals, such 
that G(xt,xY) = </> i (x)Gi i (t,^)^(x')- The Green 
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FIG. 2: The correlation part of the second-order self-energy 
(a), a diagrammatic representation of the Bethe-Salpeter 
equation (b), and the corresponding Bethe-Salpeter kernel 
(c). The Green function lines represent self-consistent, full 
Green functions. 



function is consequently represented as a time-dependent 
matrix, and the equations of motion reduce to a set of 
coupled matrix equations. The HF orbitals <pi{x) are 
themselves given by linear combinations of Slater func- 
tions centered on the nuclei of the molecule. 

We have solved the Kadanoff-Baym equations within 
the second-order self-energy approximation, illustrated 
in Fig. [2] This is a conserving approximation [T" 
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which is essential for obtaining results in agreement 
with macroscopic conservation laws. This is easily ver- 
ified numerically, for instance by checking that the to- 
tal energy remains constant when the Hamiltonian is 
time-independent. For a given G(t,t') matrix, the self- 
energy matrix (for a spin-unpolarized system) is given by 
E(t,f) = *(t,f) sHF (*) + £ (2) (*,*')> where 

= J2 G kl (t,t')G mn (t,t')G pq (t',t) 

klrnnpq 



and E5 F (*) = -iE W G kl (t,t+)(2i 



U jk ) is the HF 



self-energy. The two-electron integrals are defined by 



v m = / / rfxdx>*(x)^(x')t;(r-r')^(x')^(x). (3) 



As the initial Hamiltonian is time- independent, the 
Green function on the imaginary track of the contour 
only depends on the difference between the two imag- 
inary time-coordinates. Solving for G M (it — it') = 
G(—ir, —it 1 ) (we use r to denote time-arguments on the 
imaginary axis) is then equivalent to solving the ordinary 
Dyson equation within the finite temperature formalism 
0. Note that the definition of G M in Ref. differs 
from the one used here by a prefactor i. Since the self- 
energy is a functional of the Green function, the Dyson 
equation should be solved to self-consistency. 



Once the equilibrium Green function has been calcu- 
lated, it can be propagated according to the Kadanoff- 
Baym equations Q. Time- propagation means that the 
contour which initially goes along the imaginary axis is 
extended along the real axis (see Fig. [T]). The Green 
functions with both arguments on the real time-axis 
are denoted by the functions GfAt,t') = —i(ci(t)Cj(t')) 
(if t is later on the contour than t') and Gfj(t,t') = 
i{cAt')ci(t)) (if t' is later than t), with the symmetry 



Gfj {t, t) 



G ij 



(t,t) 



-Gfiiff ,t) and the boundary condition 



G<(t,t) 
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We also need to calculate 



the functions G^(t,ir) and G^(ir,t) with one real and 
one imaginary time- argument. The implementation of 
the propagation is similar to the scheme described by 
Kohler et. al. in Ref. [13], with one important difference 
being that we are here dealing with an inhomogenous 
system, i.e. the Green function, self-energy, and h(t) are 
time-dependent matrices rather than vectors. 

Another important difference with the propagation 
scheme described in Ref. |17| is the initial correlations. In 
our case, they are given by the equilibrium Green func- 
tion according to G<(0,0) = G M (0~) and (0, -it) = 
G M (-ir). Due to the anti-periodicity G M (ir + if3) = 
—G M (ir), the resulting nonequilibrium Green function 
will automatically satisfy the Kubo-Martin-Schwinger 
boundary condition G(0,t) = —G(—if3,t). The time- 
stepping follows the four coupled Kadanoff-Baym equa- 
tions 

[id t -h(t)]G > (t,t') = I>(t,t'), 
[id t -h{t)]G\t,iT) = lXt,ir). (4) 

and the two corresponding adjoint equations [16]. The 
collision integrals are 

I?(t,t') = f dii; R (t,t)G > (i,t')+ f diZ > (t,i)G A (i,t') 
Jo Jo 

i r 

+- / df tf(t,-i?)G l (-i?,t') (5) 
1 Jo 

I ] (t,ir) = [ dtZ R (t,t)G ] {liT) 
Jo 

1 f 8 

+ - / dfj:' ] (t,-if)G M (if -ir). 
i Jo 
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The retarded and advanced functions G R I A and Y, R / A are 
defined according to F R / A (t,t') = 5(t - t')F s {t) ± 6(t - 
t') [F> (t, t')-F< (t, t')} , where only the self-energy has a 
singular part (the HF self-energy). The last terms in each 
of Eqs. ([5|) and © account for the initial correlations of 
the system, and do not vanish when t,t' — ► 0. 

We have been able to propagate the Green function 
for a number of closed-shell atoms and small diatomic 
molecules, where one can aim at quantitative agreement 
with experimental results. The kind of calculations pre- 
sented in this paper typically take 48 hours. The basis 
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set must be large enough to describe the essential details 
of the electron dynamics. The most important limiting 
factor is the energy level structure of the systems; for 
heavier atoms, the large eigenenergies of the core levels 
lead to rapid oscillations in the Green function, and one 
consequently needs to propagate using time-steps much 
smaller than the time-scale of the interesting physical 
phenomena dominated by the valence electrons. We have 
in these calculations used 28 basis functions for beryllium 
and 25 functions for H2, with orbital energies lower than 
5 Hartree. We include a time-dependent electric field 
in the direction of the molecular axis, so that the sys- 
tem preserves a cylindrical symmetry. Generalizing this 
scheme to systems of lower symmetry does not lead to 
other complications than increasing the size of the cal- 
culations. Since the systems considered here only have 
discrete energy levels, we do not observe strong damp- 
ing effects such as what is observed in systems with a 
continuous spectrum 17, HI]. 

Figure [3] shows the imaginary part of Tr{G < (ti, £2)}, 
calculated for an H2 molecule in equilibrium, and in the 
presence of a constant electric field E(t) = 8(t)Eo di- 
rected along the molecular axis. In both cases, the trace 
of G < along the time-diagonal is constant, and equal to 
the particle number. In the ground-state (shown to the 
left), the Green function only depends on the difference 
<i — t%, and the oscillations perpendicular to the time- 
diagonal are given by the ionization potentials of the 
molecule. The right figure illustrates how the electric 
field (E = 0.14 a.u.) changes the spectral properties 
of the molecule. In addition to the expected narrow- 
ing of the ridge along the time-diagonal, the oscillations 
along the t\ — ti direction are damped. In the upper 
figures, we show ImTr{G < (ii, t 2 )} for a fixed T = (ti + 
i2)/2 = 20 a.u., compared with the same function cal- 
culated from time-dependent HF (TDHF). The ground 
state HF Green function has the form Tr{G^ F (ii, £2)} = 



and for the H2 molecule we therefore 



haveImTr{G^ F (T + t/2,T-i/2)} = 2cos(eii). The cor- 
related Green function has a sharper peak at t = and 
while it is periodic in the relative time-coordinate, it is 
characteristic of a correlated spectral function that it can 
not be fitted to a cosine function at small t [1(3] ■ With an 
added electric field, the energy fed into the system leads 
to a narrowing of the spectral function peak. 

Time-propagation is also useful as a direct method for 
calculating response functions and excitation energies [f| . 
The excitation energies of the system can be obtained 
from the poles of the density response function x{ u )i de- 
fined by Sn(u>) = x R ( LU )fi v ( UJ )- Perturbing the system 
with a "kick" of the form 5v(t) = V5(t) excites all states 
compatible with the symmetry of the perturbing poten- 
tial V. For a kick in the form of an electric field along 
the molecular axis, the induced dipole moment is given 
by d(t) — Eq J drdr' zx fl (r, r'; t)z'. The imaginary part 
of the Fourier transformed dipole moment then has peaks 
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FIG. 3: The lower figures show the trace of the imaginary 
part of G < (t\,t2) for an H2 molecule in its ground state (left) 
and in an applied electric field (right). The Green function 
in equilibrium only depends on ti — t%. The upper figures 
shows the same quantity for a fixed value of T = (t\ + tz)/2, 
compared with the same function obtained from TDHF. 




FIG. 4: The imaginary part of the polarizability a(uo) of a 
beryllium atom calculated from the Green function propa- 
gated up to a time T. 



at the poles of x fl (w), corresponding to the excitation en- 
ergies of the system. Time-propagation is in this way an 
interesting, and far more direct alternative to calculating 
the response function from x(F2) = L(l,2;l,2) where 
the particle-hole propagator L is found by solving the 
Bethe-Salpeter equation L = Lq + LqKL, as 

illustrated diagrammatically in Fig. [2] The self-energy 
approximation used in this paper would correspond to 
the kernel K = ST,/ 5G shown in Fig. [2j where it should 
be noted that the Green functions are the full Green func- 
tions of the interacting system. In Fig. 01 we have plot- 
ted the imaginary part of the polarizability, defined ac- 
cording to ar{i-o) = —1/Eq J q T dt e lujt d(t) , of a beryllium 
atom for various durations T of the time-propagation. 
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be necessary to use pseudo-potentials in order to avoid the 
very short time-steps necessary to account for the core- 
levels. For larger molecules or solids, the calculations 
are certainly feasible in the framework of model hamil- 
tonians that include electron interaction. Green function 
calculations are therefore highly important 1) for provid- 
ing benchmark results for simpler methods, 2) for inves- 
tigating the role of electron correlation, and 3) as and 
alternative to solving the Bethe-Salpeter equation with 
highly sophisticated kernels. 



FIG. 5: The time-dependent dipole moment of a Be atom, 
calculated within the HF and the second-order approxima- 
tion. 



The polarizability develops a distinct peak at the S — ► 
1 P transition energy lo = 0.189 a.u. (compared to the 
TDHF value of 0.178 calculated within the same basis, 
and the experimental value of 0.194 [l9(), which becomes 
increasingly sharp as the propagation time is extended. 
As the system consists only of discrete energy levels, the 
damping of the time-dependent dipole moment d(t) as a 
function of time is not significant in the relatively short 
duration of the time-propagation. The position of the ex- 
citation energy peak in Fig. [4] therefore converges slowly, 
but extrapolation schemes nevertheless gives the position 
of the peak very accurately. 

In Fig. [51 we have plotted the time-dependent dipole 
moment of the beryllium atom, this time with a non- 
perturbative "kick" at Eq = 1.0 a.u.. In this case, we 
see a clear difference between the uncorrelated Hartrce- 
Fock result and the dipole moment calculated from the 
Kadanoff-Baym equations. The difference between the 
correlated and the uncorrelated results become more ap- 
parent with time, due to the non-linear effects entering 
via the self-energy. Since we have now excited all unoc- 
cupied states, we can expect some damping in the time- 
dependent dipole moment, but the damping is much more 
pronounced in the correlated calculation. 

In conclusion, we have demonstrated that time- 
propagation of the non-equilibrium Green function can 
be used practical method for calculating non- 

equilibrium and equilibrium properties of atoms and 
molecules (e. g. for atoms in strong laser fields). For 
these small systems, the second order approximation is 
clearly well-suited. For extended systems, such as molec- 
ular chains, it becomes essential to cut down the long 
range of the Coulomb interaction. This is done effec- 
tively within the GW approximation [2^, [2l[ (known as 



the shielded potential approximation in Rcfs. [13, |l4j) 
where the self-energy is given as a product of the Green 
function G and the dynamically screened interaction W 
and which we have currently implemented for the ground 
state [22[ . For calculations involving heavier atoms it will 



[1] Y. Xue, S. Datta, and M. A. Ratner, Chem. Phys. 281, 
151 (2002). 

[2] L. V. Keldysh, Zh. Eksp. Teor. Fiz. 47, 1515 (1964), [Sov. 
Phys. JETP, 20, 1018 (1965)]. 

[3] L. P. Kadanoff and G. Baym, Quantum Statistical Me- 
chanics (W. A. Benjamin, Inc., New York, 1962). 

[4] W. Schafer, J. Opt. Soc. Am. B 13, 1291 (1996); 
M. Bonitz, D. Kremp, D. C. Scott, R. Binder, W. K. 
Kraeft, and H. S. Kohler, J. Phys. Condens. Matter 8, 
6057 (1996); R. Binder, H. S. Kohler, M. Bonitz, and 
N. Kwong, Phys. Rev. B 55, 5110 (1997). 

[5] N.-H. Kwong and M. Bonitz, Phys. Rev. Lett. 84, 1768 
(2000). 

[6] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 
(1984). 

[7] M. A. L. Marques, C. A. Ullrich, F. Nogueira, A. Rubio, 
K. Burke, and E. K. U. Gross, Time- dependent density 
functional theory (Springer, Berlin Heidelberg, 2006). 
[8] R. van Leeuwen, Phys. Rev. Lett. 76, 3610 (1996). 
[9] H. Haug and A.-P. Jauho, Quantum kinetics in transport 
and optics of semiconductors (Springer, Berlin Heidel- 
berg, 1998); N.-H. Kwong, M. Bonitz, R. Binder, and 
H. S. Kohler, Phys. Stat. Sol. 206, 197 (1998). 

[10] M. Bonitz, Quantum Kinetic Theory (B. G. Teubner, 
Stuttgart; Leipzig, 1998). 

[11] P. Danielewicz, Ann. Phys. (N. Y.) 152, 239 (1984). 

[12] We here consider only finite systems, for which the finite- 
temperature formalism is used as a formal device for im- 
posing boundary conditions on the Green function, but 
is only meaningful in the T — > limit. 

[13] G. Baym and L. P. Kadanoff, Phys. Rev. 124, 287 (1961). 

[14] G. Baym, Phys. Rev. 127, 1391 (1962). 

[15] N. E. Dahlen and R. van Leeuwen, J. Chem. Phys. 122, 
164102 (2005). 

[16] N. E. Dahlen R. van Leeuwen and A. Stan, J. Phys.: 

Conf. Ser. 35, 340 (2006). 
[17] H. S. Kohler, N. H. Kwong, and H. A. Yousif, Comp. 

Phys. Comm. 123, 123 (1999). 
[18] D. Semkat, D. Kremp, and M. Bonitz, Phys. Rev. E 59, 

1557 (1999). 

[19] A. Kramida and W. C. Martin, J. Phys. Chem. Ref. Data 

26, 1185 (1997). 
[20] L. Hedin, Phys. Rev. 139, A796 (1965). 
[21] F. Aryasetiawan and O. Gunnarsson, Rep. Prog. Phys 

61, 237 (1998). 

[22] A. Stan N. E. Dahlen and R. van Leeuwen, Europhys. 
Lett. 76, 298 (2006). 



